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Abstract 

A non-reflecting boundary condition (NRBC) 
for practical computations in fluid dynamics 
and aeroacoustics is presented. The technique 
is based on the hyperbolicity of the Euler equa- 
tion system and the first principle of plane 
(simple) wave propagation [1]. The NRBC is 
simple and effective, provided the numerical 
scheme maintains locally a C 1 continuous so- 
lution at the boundary. Several numerical ex- 
amples in ID, 2D and 3D space are illustrated 
to demonstrate its robustness in practical com- 
putations. 

1 Introduction 

It is well-known that non-reflecting boundary con- 
ditions (NRBCs) play an important role in fluid flow 
and aeroacoustics computations. The need for artificial 
boundary conditions arises when the domain of the prob- 
lem is unbounded and extends to infinity. In order to 
treat the problem numerically, a domain of finite size is 
required and artificial boundaries are imposed. At these 
artificial boundaries, NRBCs are sought for to minimize 
their influences on the flow. A spurious reflection result- 
ing from an inappropriate numerical boundary condition 
will contaminate the flow field and may entirely spoil 
the flow computation. Research on NRBC is a challeng- 
ing topic in engineering and applied mathematics. For 
decades, a vast number of papers on NRBC have been 
published, e.g. see [3-6], the review paper by Givoli [2], 
and the references cited there. 

In one-dimensional flow, at an artificial boundary, En- 
quist and Majda [3], and Hedstrom [7] proved that a 
boundary condition is non-reflecting is equivalent to say- 
ing that the characteristic variables corresponding to the 
incoming characteristic curve remain constant across the 
artificial boundary (see also Hirsch [8], p. 370). For 
multi-dimensional flow, this 1-D technique is combined 
with dimension splitting and applied in the practical 


NRBC treatment. Such a combined treatment has been 
the topic of many papers on the characteristic based NR- 
BCs (see e.g. [4] and references in [2]). 

Other alternative treatments for NRBC found in the 
literature include (see [2-6]): 

(i) in the far-field, using a predictable asymptotic ana- 
lytical solution at the boundary ( the radiation boundary 
condition), 

(ii) diminishing the strength of the waves/disturbances 
before they reach the artificial boundary, and thus mini- 
mizing the reflecting effect. Usually, increased numeri- 
cal damping is applied to a zone between the core domain 
and the artificial boundary (the buffer zone or sponge 
zone) to do the job. In the recently developed PML (per- 
fectly matched layer) method, a specially designed equa- 
tion system is imposed in the matching layer (or sponge 
zone) to guarantee the exponential decaying of the dis- 
turbances in the layer [5,6]. 

In the present paper, a different but simple criterion 
is introduced to treat the NRBCs of the time-dependent 
hyperbolic conservation laws of gas dynamics. The cri- 
terion is based on the first principle of plane (simple) 
wave propagation [ 1 ] rather than the characteristics the- 
ory. Emphasis is put on their viable practical applica- 
tions. As it turns out that the NRBCs used in the recent 
CE/SE finite volume schemes for flow and aeroacous- 
tics computations (e.g. [ 14,15]) can be directly derived 
from this criterion, the present paper also serves to ex- 
plain why these simple NRBCs of the CE/SE schemes 
work. 

The paper is arranged as follows: In Section 2, based 
on the hyperbolicity of the Euler equation system and 
the propagation of plane wave, the continuity criterion 
of NRBC is introduced and proved. An extrapolation- 
like NRBC (Type I) based on this criterion and the nu- 
merical procedure are described. Then the relation be- 
tween the NRBC and the flux balance across the bound- 
ary surface is established, which leads to another type 
of NRBC (Type II). In Section 3, several numerical ex- 


1 

American Institute of Aeronautics and Astronautics 


amples for outflow NRBC in one and multi-dimensional 
space are presented. Numerical examples with Type II 
NRBC at the inflow and other artificial boundaries are 
demonstrated in Section 4. Application of buffer/sponge 
zones is illustrated in Section 5. At last, the paper is con- 
cluded with remarks in Section 6. 

As the time-dependent hyperbolic conservation laws 
of gas dynamics ( in dimensionless form) is always in- 
corporated in the new treatments of NRBC in the present 
paper, they are briefly described here for later use: 

Ut 4- F x + G y + H z = Q, (1) 

where x, y , 2 and t are the streamwise and transversal co- 
ordinates and time, respectively. The conservative flow 
variable vector U and the flux vectors in the streamwise 
and radial directions, F, G , and H , are given by: 
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where u , v, w and p, p are respectively the velocity com- 
ponents, density and pressure, e = p (y-i) + l/2(^ 2 + 
v 2 -I- w 2 ), and the enthalpy H = p/p + e with 7 = 1.4. 
The right hand side Q is the source term which may in- 
clude the possible external forcing terms and/or viscous 
fluxes. 

By considering (x,y,z,t) as coordinates of a four- 
dimensional Euclidean space, £4, and using Gauss’s di- 
vergence theorem, it follows that Eq. (1) is equivalent to 
the following integral conservation laws: 



where S(V) denotes the surface around a volume V in 
£4 and I m = ( F m ,Gm,H m ,U m ) stands for the flux 
vectors, I m ds = I m • n ds , n being the outgoing unit 
normal vector in £4, is the flux at the infinitesimal sur- 
face element ds. 



2 The continuity criterion of NRBC for time 
dependent conservation laws 

There are various techniques to treat the NRBC based 
on characteristics theory of the hyperbolic system. For 
instance, a well-known 1 -D flow NRBC treatment by En- 
quist and Majda [3], and Hedstrom [7] is the requirement 
that the local perturbation (disturbance) along incoming 
characteristics be made vanish at the boundary (see [8], 
p.370). Let W = ( W\,W2,W3) t be the 1-D character- 
istic flow variables, the above requirement states that for 
those k such that the corresponding characteristic enters 
the computational domain through the artificial bound- 
ary: 

A w k = 0 (3) 

Thompson [4], for example, has details on the practical 
implementation of the characteristic NRBCs with a finite 
difference scheme. Hirsch [8] also offers an excellent 
resource of various NRBCs. 

In finite volume numerical approaches with hyperbolic 
conservation laws, grid nodes are often cell centers and 
the boundary faces are often formed by the boundary 
cell surfaces. No node lies exactly on the boundary. As 
such, a continuity criterion of NRBC and the consequent 
NRBC treatments are recommended based on the hy- 
perbolicity of the equation system and the first principle 
of plane wave propagation [ 1 ]. They are simple, robust 
and particularly appropriate for cell center finite volume 
schemes. Their limitations are also briefly discussed. 

2.1 The continuity criterion of NRBC 

The continuity criterion is first introduced in an heuristic 
way and then proved via the first principle of plane wave 
propagation [1]. 

Under a general assumption that the flow is continu- 
ous near the boundary, i.e., with no shock or contact dis- 
continuity, we first consider the behavior of the charac- 
teristic variables W = (wuW 2 ,Ws,W 4 ,Ws) T across the 
artificial boundary surface element As as time elapses. 
As also represents the interface between a boundary cell 
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and its corresponding ghost cell. Note that any (artifi- 
cial) spatial boundary is a cylindrical hyper-surface in the 
space-time £4 space (Fig. 1 demonstrates the situation in 
£3). 

Let Pi and P e be respectively an interior point and 
an exterior point of the computational domain in the £4 
space. Both Pi and P e lie in the neighborhood of (9, a 
point on As (Fig. 1). Recall the standard NRBC treat- 
ment for 1-D flow [3,7], 


Aw k = w k (Pi) — Wk{P e ) = 0 , for selected fc. 




Figure 2: Sketch of the continuity criterion in 1 -D flow, 
only a real component of V is shown. 


When Pi tends to O from the interior of the domain and 
P e tends to O from the exterior of the domain, the NRBC 
Eq. (3) becomes: 

lim w fc (Pj) = Jim u/k(P e ) = w k (0) (4) 

Pi-tO Pe — tO 

for those selected fcs. Therefore, the usual NRBC treat- 
ment is formally interpreted as a continuity problem of 
w k across the boundary surface. 

Secondly, for the NRBC of multi-dimensional flow, 
we formally extend the continuity of w k across the 
boundary surface to all fcs rather than just a selected num- 
ber of fcs: 

lim Wk(Pi) = lim w k (P e )> for all k. (5) 

Pi -yO P e —yO 

Let U = (ui,u 2 ,U 3 , u 4 ,us) t = {p,pu,pv,pw,pe) T 
and V = (vi,V 2 ,V 3 ,V 4 ,v 6 ) t = ( p,u,v,w,p) T here- 
spectively the conservative flow variables and the prim- 
itive flow variables. From the local equivalence of 
the characteristic variables W and V or U (see e.g. 
Hirsch[ 8 ], p. 1 55- 1 56), Eq. (5) is equivalent to the conti- 
nuity of V or U across the artificial boundary surface: 

lim u k {Pi) = lim Uk(P e ), for all fc ( 6 ) 

Pi-yO Pe-yO 

or 

lim Vk(Pi) = lim v k (P e ), for all fc (7) 

It is advantageous to switch from the continuity relation 
(5) of W to that of V or U, ( 6 ) or (7), since the latter 
can be treated in an easy way. At this stage, the conti- 
nuity criterion of NRBC is heuristically inferred and will 
be proved later. The feasibility that all u k s (v k s and w k s 
as well) can be made approximately continuous simulta- 
neously at the boundary surface is demonstrated in the 
numerical examples of NRBC in § 2.2. 1 . 

Now, with the presence of the artificial boundary s 
(hyper-surface), £4 is bisected into two portions, domain 
interior Di and domain exterior D e . Within each portion, 
the flow is governed by the same Eq. ( 1 ). From the first 
principle of plane wave propagation, it can be shown that 


locally there is no reflection at the artificial boundary sur- 
face s if the continuity criterion Eq. (7) (or Eq. ( 6 ), or 
Eq. (5) ) is satisfied. 

Proof: Consider a non-conservation form of Eq. ( 1 ) 
in primitive variables V: 


dV ~<9V 
dt + dx 


ay oz 




( 8 ) 


where A , £, and C are the jacobian matrices and func- 
tions of V. Q is the source term vector. As a result of 
the continuity condition (7), an admissible given set of 
V at the boundary s may be used as a common boundary 
condition to solve separately for V 7 and V e in their cor- 
responding subdomains Di and D e , which are separated 
by the artificial boundary s. (Note that generally, the ad- 
missible V given at s should be identical to the solution 
of V over the entire domain, see Appendix II). Here, V 7 
and V e are respectively the solutions of ( 8 ) in Di and 
D e . Let V be the solution of ( 8 ) over the entire domain. 
Due to the uniqueness of solution for well-posed initial- 
boundary value problems (Appendix II), V; is identical 
to V in Di and V e is identical to V in £ e . Therefore, 
in a neighborhood of s, V* and V e are mutually a con- 
tinuation of each other across the boundary s and hence 
there is no reflection (Fig. 2). 

To be more specific in terms of plane wave propaga- 
tion, notice that from the hyperbolicity of the system ( 8 ), 
for any propagation direction k = (fc x , fc y , k z ) (the wave 
number vector) , the matrix K — k x A + fc y B\-k z C has 
real eigenvalues and linearly independent left eigenvec- 
tors, and then, there exists a plane wave (or simple wave) 
solution: 

V — y e *(k«x-w t) (9) 

where x =(x,y,z) is the position vector, i = >/“X w 
is the angular frequency that is related to the eigenval- 
ues via the dispersion relation (see Courant and Hilbert 
[1], and Hirsch [ 8 ], p. 147, 150, and Appendix I). As a 
wave solution of the hyperbolic system ( 8 ) can be locally 
written as a superposition of the plane wave solutions by 
Fourier integral (Appendix I), it suffices to consider only 
the behavior of a single plane wave solution in the form 
of (9) at the artificial boundary s. 
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boundary 



Figure 3: Numerical treatment of NRBC in E 3 . 

Let O (x 0 ,t 0 ) be any point at the artificial boundary 
s, then ( 8 ) can be locally linearized in the neighborhood 
of O, i.e., A, B and C are frozen at V Q . For any given 
wave number vector k, from the continuity criterion (7), 
V 0 = Vi = V e , all A, B and C remain the same across 
the boundary and so also the eigenvalues cjj = w e = u) 
(see Appendix I). At O, the plane waves V; and V e share 
the same k,x, a ; and t, and hence the same phase. Again, 
from the continuity criterion (7), V* = V e , or 

_i(k«x-w t) i(k*x-a> t) 

,c — V e t, , 

i.e., Vj = V e . Thus, the plane waves V* and V e share 
the same amplitudes too. Therefore, V* is completely 
identical to V e in a neighborhood of O at the boundary 
(in terms of phase and amplitude), and there is no reflec- 
tion at O across the artificial boundary surface. 

The continuity of V, or U or W (Eq.(5)- Eq.(7)) 
across the boundary surface is thus the basic criterion of 
NRBC adopted in the present paper. In §2.2, the numeri- 
cal NRBC (Type I) is constructed based on the continuity 
criterion. In §2.3, the relation between an NRBC and the 
flux balance across the boundary surface is established. 
Such relation leads to another absorbing NRBC (Type 

II). 

2.2 The numerical treatment of NRBC 

Fig. 3 illustrates a 2-D (in E 3 ) NRBC treatment. Let 
A ABC be a boundary cell centered at P, with the 
side BC coincident with the artificial boundary surface. 
A BCD is the ghost cell centered at Q , sharing the 
boundary edge BC with A ABC. Let O be the centroid 
of the boundary surface element BCC' B' . The limiting 
process of limp._>o U (Pi) is equivalent to extrapolating 
U from the interior node P to O by Taylor expansion. 
Similarly, limp e _>o U(P e ) is equivalent to extrapolating 
U from the exterior ghost node Q to O by Taylor expan- 
sion. 


Although theoretically, ( 6 ) implies up to C°° continu- 
ity, in numerical approximation, only low order continu- 
ity such as C°, C l or C 2 , etc. can be achieved. Since 
a plane wave solution (9) is based on two parameters, its 
amplitude and phase, the numerical approximation of U 
(or V) is required to be at least C 1 continuous at the arti- 
ficial boundary in order to be consistent with the physical 
solution. Taking a 1 -D version of (9) for example, the C l 
continuity requirement is explained as follows. 

It suffices to consider only a scalar component of V, 
say, the first one p. After discretization, the (artificial ) 
boundary surface element center O is used to represent 
the entire surface element As. Then, from the continuity 
criterion, approximately, it can be inferred that the ‘am- 
plitude’ ( p 0 )i = ( p 0 )e — Po, where the subscripts 0 , 
i and e stand respectively for the surface center O, do- 
main interior and exterior. Approximately, at the bound- 
ary surface As, 

Pi = p 0 e i( - kiXo ~ Wit °) = p e = p 0 e i( - k ‘ Xo ~ u ‘ to K ( 10 ) 

Note that numerically the C° continuity result (10) pro- 
vides no information about the wave number k and the 
frequency cj. With the presence of phase error, numeri- 
cal reflection may still occur. However, if the numerical 
continuity is enhanced from C° to C 1 , i.e. 

( Pi)x = ikipi = ( Pe)x = i^ePey 
( Pi)t ~ i^ipi — ( Pe)t — 

then ki = k € and = c u e and there is no phase error. 

In constructing the NRBCs, although any one of U, V, 
and W can be used, U is selected since it is employed in 
the numerical examples in the present paper. Therefore, 
in addition to U p , the space and time gradients of U at 
P , namely, U x , U y , U z and U t are also required. The 
resulting linear Taylor expansion (C 1 continuity) yields 
better accuracy and is consistent with the NRBC crite- 
rion. The NRBC at the ghost node Q now turns out to 
be a problem of how to define U and its gradients at Q 
so that the flow is C 1 continuous at the boundary surface 
(represented by O). 

2.2.1 Examples of NRBC For the Type I (out- 
flow) NRBC, under a mirror image assumption explained 
later, it is found that a simple extrapolation technique 
works well. 

First, an example of NRBC in E 3 (2-D space) for tri- 
angular mesh is illustrated. As shown in Fig. 3, assume 
A ABC is a triangular boundary cell with the edge BC 
lying on the boundary and conveniently parallel to the y- 
axis. Define a ghost node D as the mirror image of the 
triangle vertex A with respect to BC. Then A ABC and 
A DBC are mutually mirror images of each other (the 
mirror image assumption). At time step n, conservative 
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variables U are given at the cell center P of A ABC. 
Then, the NRBC (labeled as Type I) at the geometrical 
center Q of the ghost cell may be defined as: 


(U)q = (U)p, (U X ) Q = (U x )p = 0, (U y ) Q = (U y )p. 

( 11 ) 

Apply linear Taylor expansions to domain interior and 
exterior separately: 


(U o)interior = Up + {VO “ 2/p)(U y )p + 1/2 At(U t ) p. 



(U o)exterior — Uq -f- (t/o 2/Q)(U y )q ~f* 1/2 A^(U(;)q 

/iCTlce, (U o)exterior = {VJ o) interior • 

Here, the subscripts P, Q and O of x and y denote the 
corresponding coordinates of P, Q and 0, and from ( 1 ), 
the time partial derivatives (Ut)g = (U t )p can be di- 
rectly obtained. Thus, U is C l continuous at O across 
the boundary surface element, (6) is satisfied and the 
boundary surface element is non-reflective. 

In a consistent way, for 3-D flows, under the same mir- 
ror image assumption on the ghost cells, the following 
extrapolations are valid NRBCs with C l continuity: 

U q = U p ,(U x ) q = (U x )p-0, 

(U y ) Q = (U y )p,(U z ) Q = (U z )p, (12) 

or 

U Q - Up + Ax(U x )p, (U X ) Q = (U X ) P , 

(U y ) Q = (U y )p,(U z ) Q = (U z ) P , (13) 

where Ax = xq — xp. 

As demonstrated in the examples in §3 and §4, this 
Type I NRBC works well for either supersonic or sub- 
sonic flows at the outflow boundary, but it should be 
noted that: 

1 . ( 1 2) or ( 1 3) is but a possible selection under the mir- 
ror image assumptions, there are other forms of NR- 
BCs based on (5) - (7); 

2. The extrapolation technique utilizes the nearby U p 
data to approximate the U data at the artificial 
boundary, which is not an unreasonable choice, but 
there is a danger that the solution could drift away 
from the true solution (see Appendix II). A remedy 
is to incorporate the Type I NRBC with other phys- 
ical boundary conditions (e.g. back pressure, etc.) 


Figure 4: Control volumes ( CVs ) in P 3 for compact 
updating. Left: boundary cell A ABC and the corre- 
sponding hexagon CV base ASBQCRA ; right: quadri- 
lateral boundary cell ADCD and its corresponding oc- 
tagon CV base ATBQCRDSA ; R , 5, T are centers of 
neighboring cells and BC the boundary. In any case, 
quadrilateral PBQC is a portion of the CV base. 

2.2.2 Numerical implementation of NRBC 

The implementation of the NRBC is incorporated in the 
numerical procedure and may be summarized in the fol- 
lowing steps: 

(i) Based on the flow data at the boundary cell center 
P, i.e. Up and its spatial gradients (slopes) U x , U y and 
U z , determine the flow data Uq as well as its gradients 
U x , U y and U z , i.e. the NRBC at Q , as described in 
§ 2 . 2 . 1 . 

(ii) Update U at boundary cell center P to the new 
time level by the conservation laws (2). In order that U 
be C 1 continuous across the artificial boundary, the up- 
dating procedure must be carefully designed to take ac- 
count of the accuracy of surface flux calculation. Here 
a compact updating procedure described in §2.2.3 is rec- 
ommended. 

(iii) After U at all the interior cell centers of the com- 
putational domain are updated, evaluate the new spatial 
gradients U x , U y , U z at the boundary cell center P by 
finite difference. For multi-dimensional flows, a linear 
equation system is required to solve for the gradients. 

(iv) Repeat steps (i) - (iii) and march in time. 

2.2.3 Compact updating The updating proce- 
dure recommended here is identical to the one used in 
the recent CE/SE method (Chang et a/ [10,1 1 ]) and sim- 
ilar to the NT (Nessyahu-Tadmor ) scheme [12] in 1-D 
flow. 

The purpose for compact updating is to achieve high 
accuracy (C 1 continuity) from a small cell stencil (e.g. 
the stencil formed by the immediate neighboring cells). 
Therefore, not only U but also its gradients U x ,U y 
and U z are required. The compact updating is capable 
of maintaining C 1 continuity of U across the artificial 
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control volume 



Figure 5 : Fluxes balance on an internal face Si in £ 3 . 


boundary. 

Consider a triangular cell A ABC in £ 3 (Fig. 4 , shaded 
area). Let P and Q,R,S be respectively the cell centers 
of AABC and its neighboring triangular cells. Most of 
the finite volume schemes use the space-time cylinder 
AABC as the control volume ( CV ) for updating V or 
U at P. Surface flux along, say, AB is obtained by ex- 
trapolation from 5 across the cell to the surface center 
along AB. Another treatment is to include S in the CV 
and replace the flux along AB by the fluxes along ,45 
and SB. Since 5 sits right on both surfaces along AS 
and SB , no extrapolation across the interior of the CV 
is involved. For the triangular cell A ABC in Fig. 4 the 
CV turns out to be a space-time hexagon cylinder in £3 
based on ASBQCRA. By applying the integral conser- 
vation laws (2) to the CV, updating flow data at P based 
on the flow data at a compact node stencil of Q,R,S 
is now completed. Fig. 4 also demonstrates that for 
quadrilateral mesh cells the CV is an octagon cylinder 
in £3 ( 2 -D space). 

The surface fluxes for the CV surfaces passing through 
a vertex, say, cell centers Q, R, 5 , can be evaluated 
by first extrapolating U along the surface to their cor- 
responding surface centers by linear Taylor expansion, 
calculating flux functions F Q and H, and then incor- 
porating the surface unit normal vector and computing 
the fluxes. For high space dimensions, the CV s are geo- 
metrically more complicated. More details including the 
updating of U x , U y , and U c can be found in [10,1 1]. 

It should be noted that the compact updating is sug- 
gested for boundary cells only. For interior cells, one 
may still continue to use other finite volume schemes. 

2.3 Relation between an NRBC and the flux 
balance across the boundary surface 

In this subsection, the relation between two statements is 
established. The first one states that the incoming fluxes 
at the artificial boundary surface are equal to the outgo- 
ing fluxes, or fluxes are balanced across the boundary 


Figure 6: Fluxes balance across a boundary surface ele- 
ment in a control volume (the above figure) in £3. 

surface. Here, the outgoing flux is defined as a portion of 
the left hand side of (2): 



where A 5 denotes the artificial boundary surface (el- 
ement). The second statement is that the boundary is 
non-refiective. Such a relation leads to the Type II 
NRBC. Numerical implementation of the Type II NRBC 
is straightforward and will be illustrated in § 2 . 4 . 

Let V be any control volume in the £4 space inter- 
sected and divided into two portions V\ and V2 by an 
internal surface Si . Let £1, £2 and < 7 i, <72 be the fluxes 
around the surfaces of V\ and V2 respectively. Here <j\ 
and <72 are the outgoing fluxes at the interface Si for V\ 
and V2 respectively (Fig. 5 ). Then, the following lemma 
holds: 

Lemma 1 : For a control volume V in the £4 space, 
fluxes passing through any of its internal surface 5 * are 
balanced, he. <71 + 02 = 0 . 

Proof: Apply ( 2 ) to V, V\ and V2 separately. We have 

£1 -b £2 = / QmdV, 

Jv 

£l+< 7 i= / QmdV, £2 + <72 = / Q m dV. 

J V\ Jv 2 

Add the second and the third equations and then subtract 
the first one, we have g\ + <72 = 0 . 

In other words, across any interior surfaces, no extra 
‘source’ should be generated. 

Next, we show that locally the continuity of flow vari- 
ables (Eq. ( 5 ), ( 6 ) or ( 7 )) can be inferred from flux bal- 
ance. Hence, flux balance suffices for an NRBC. 



artificial boundary surface 
element ds 


control volume 
V=Vi +V 2 
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Consider an element ds of the cylindrical spatial 
boundary surface in E4. As shown in the control vol- 
ume in Fig. 6 (in E$), ds is centered at O. As- 
sume the outgoing unit normal vector at O is ni = 
(n x ,n y ,n z ,0) T Ahen, the incoming unit normal is n 2 = 
(-n X) — n y , — n z ,0) T . The outgoing flux <n for the in- 
ner boundary surface: 


Oi - 

/C7n\ 

<712 

<713 

— (i<s j 

l Fl \ 

f 2 

f 3 

+n y 

(Ci\ 

g 2 

G 3 

+n 2 

( Hl \ 

h 2 

H 3 


<7i4 

VCT15 / 


f 4 

\fJ 


G 4 

\gJ 


h 4 

\hJ 


— ds[n x F + n^G -f- n z il\. 

Let the outgoing boundary surface flux vector be: 

L — [tijF H- TiyG 4- 72 3 H], 

then 

L = ai/ds . 

Let V = (p, u,t;, w,p) T be the primitive flow variable 
vector. L may be considered as a non-linear vector func- 
tion of V. The jacobian matrix has the eigen values 
(e.g. see Hirsch [8], p. 1 77): 

Ai — A 2 — A3 — 72 * tlj "f tl ' 72 y 4" 71? * 72 z , 

A4 — Ai — c, A5 = Ai -h c. 

where c is the speed of sound. If none of the eigenvalues 
vanishes at O, the jacobian is non-singular, and an 
inverse vector function, the primitive flow variables V as 
a vector function of the surface flux L exists. Thus, 
Lemma 2: If the jacobian ^ is non-singular, then, 
locally, the primitive variables V are uniquely defined 
by the flux vector L at the center O of ds. 

Proof: Assume L(Vi) = Li and there is another V2 
in the neighborhood of Vi, such that L(V 2 ) = Li. By 
a linear Taylor expansion, 

L(V 2 )-L(V 1 ) = 

(|^)v=Vx(V 2 - Vi) + 0(|V 2 - V1I 2 ) = 0 

Hence V 2 — Vi = OorV 2 = Vi since the first order 
term cannot cancel with the second order term. 

From Lemma 2 and the NRBC continuity criterion (7), 
it is inferred that locally, under the condition that the ja- 
cobian is non-singular, the following lemma holds: 
Lemma 3: For hyperbolic conservation laws of gas 
dynamics, an element of the artificial boundary surface is 
non-reflective if its outgoing fluxes and incoming fluxes 
are equal (balanced). 


boundary surface 



Figure 7: NRBC at inflow in £3. Inflow data are given 
at the ghost cell center. 

Lemma 3 states that conditionally the flux balance 
across a boundary surface element is a sufficient condi- 
tion for NRBC. With the concept of flux balance, there 
is also an intuitive interpretation of the commonly used 
terms for NRBC such as ‘transparent’ and ‘absorbing’. 
If the flux 100% passes through the boundary surface 
element (or the flux is balanced), the boundary surface 
is said to be ‘transparent’ or ‘absorbing’ to the fluxes. 
Chang et al [9] were the first attempting to explain the 
1-D NRBC using flux concepts. 

2.4 Implementation of absorbing boundary 
condition 

The lemmas in §2.3 can be easily applied to construct 
the Type II NRBC. In case that at the ghost cell center 
nodes flow variables must be specified as the given val- 
ues (e.g. at the inflow boundary, Fig. 7), another type 
(Type II) of NRBC - absorbing NRBC arises. A control 
volume V across the boundary surface is needed to apply 
the divergence theorem (2). As shown in Fig. 7, the ghost 
cell center lies outside of the domain, and the boundary 
surface element As is an internal surface of the control 
volume V, then, from Lemma 1 and Lemma 3, As is au- 
tomatically a non-reflective boundary surface element. 

The Type II NRBC states that no extra condition of 
NRBC is needed with the prescribed flow boundary con- 
ditions at the ghost cell center. This Type II NRBC is 
flexible and valid with practically any cell shapes or con- 
figurations. Ghost cells are not required to be mirror 
images of the corresponding boundary cells. However, 
it is noted that U is still required to be C 1 continuous 
across the artificial boundary surface. Consequently, the 
compact updating is still recommended for the Type II 
NRBC. 

In many cases, over a long time period, the (subsonic) 
flow may develop a flow near the inflow boundary that is 
different from the imposed boundary conditions, or the 
boundary condition becomes ‘overdetermined’, then, a 
‘matched layer’ may be formed between the computed 
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interior flow and the imposed flow boundary condition 
(see Fig. 18). The situation is somewhat similar to 
that in the PML (perfectly matched layer) method [5,6]. 
The shock-capturing numerical scheme should be able to 
quickly resolve this non-physical discontinuity in a few 
cells. 

2.5 Discussions on the NRBCs 

In practical applications, due to discretization, approxi- 
mation and lack of information in the domain exterior, 
there are limitations for both Type I and Type II NRBCs. 

For Type I NRBC, as mentioned in §2.2.1, extrapola- 
tion technique might lead to drifting or deviation from 
the true solution because no flow data outside the out- 
flow boundary is available. In addition, under the mirror 
image assumption on the ghost cells in §2.2. 1 and the as- 
sumption that the boundary surface As is normal to the 
x axis, the NRBC (12) implies that the NRBC continuity 
criterion (5 - 7) are satisfied at any point on As. How- 
ever, as explained in the following, due to discretization 
and the possible consequent phase error, the accuracy of 
the NRBC could be degraded. 

Consider a Fourier mode in the plane wave solution 
(9): e l ( k#x_u; l \ Here, 0(x,t) = k • x — cot is the phase 
of the wave mode, with k being the wave number vec- 
tor in the propagation direction. Generally, the direction 
of k may or may not be the same as the flow direction. 
0(x, t) = const, stands for a wavefront (or a charac- 
teristic surface, see e.g. Hirsch [8], p. 150, Courant and 
Hilbert [1] ). After discretization, the center O of As 
(Fig. 3) is employed to represent the entire As. Then 
how much phase error is introduced to the Fourier mode 
by the discretization? Let x = ( x , y, z ) be the position 
vector of any point on As and xo the position vector 
of the center O. For clarity, assume time t is held un- 
changed. Then, after discretization, the phase error AO 
due to replacing x by xo is: 

AO = k*x — k*xo = k • (x — xo) (14) 

note that Ax = x — xo lies on As , AO = 0 when k 
is normal to As. Therefore, for Type I NRBC, the best 
result is obtained when the wave propagation direction is 
normal or only slightly oblique to the boundary surface. 
Otherwise, a phase error of order O(Ax) may be intro- 
duced. It deteriorates the accuracy of NRBC and causes 
numerical reflection. 

For Type II NRBC, in addition to the similar phase er- 
ror of Type I NRBC, there are other restrictions too. In 
§2.3, Lemma 3 is conditionally valid because it is based 
on the one to one correspondence between the boundary 
surface flux vector L and the primitive flow variable vec- 
tor V. The latter relies on the non-singularity of the jaco- 
bian |^. In case that is singular, Lemma 3 may fail. 
In addition. Lemma 3 is valid only locally. Globally, the 
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Figure 8: Comparison of exact and numerical results at 
the domain boundary for the 1-D Gaussian pulse prob- 
lem. 


relation between L and V involves a quadratic equation, 
the vector function V(L) could be multi-valued. This 
could break the one to one correspondence between L 
and V globally and lead to the failure of Lemma 3 in the 
global sense. 

In a nutshell, in reality, there are various situations 
that the Type I or II NRBCs can only be partially or ap- 
proximately implemented, causing spurious reflections 
at the (artificial) boundary or deviation from the true so- 
lution. In practice, an effective remedy is to impose a 
buffer/sponge zone between the boundary and the inte- 
rior domain. Although the same governing equations (1) 
or (2) are employed in the sponge zone, numerical damp- 
ing is highly increased to diminish the wave/disturbance 
amplitude before it reaches the boundary and to mini- 
mize the spurious reflection. An example is depicted in 
§ 5 . 

3 Numerical examples for outflow NRBC 

In this Section and §4, the effectiveness of these NR- 
BCs is demonstrated in numerical examples in one, two 
and three dimensional spaces. 

In principle, any finite volume scheme can be used 
with the above NRBC if it can be manipulated at cer- 
tain high accuracy. Here the recently developed space- 
time conservation element and solution element (CE/SE) 
method [10,1 1] is chosen for computing the examples 
since the compact updating is a standard procedure in 
the scheme, making application of the NRBCs straight- 
forward and effective. Full details of the method are de- 
scribed in [10,1 1]. The Type I NRBCs used with the 
CE/SE method are identical to (12), with possible minor 
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Figure 9: Propagation of 1-D Gaussian pulse at the do- 
main boundary. 



Figure 10: Free shear layer instability problem, u\ = 
Ui = 1, vi = 0, pi = 1/3.15, pi = 1, Mi = 
1.5, u 2 = U 2 = .7391 304, u 2 = 0, p 2 = 1/3.15, p 2 = 
0.5405405, with subscripts 1,2 denoting the fast and 
slow streams respectively. 


modification according to the grid layout. 


3.1 Propagation of a 1-D Gaussian pulse 

Consider a scalar initial value problem: 


du 9u_ 
dt dx 


over the range —20 ^ x ^ 450, with a Gaussian pulse 
u = 0.5exp [— (/n2)(|) 2 ] at t = 0. This is one of 
the benchmark problems of the 1st CAA Workshop [13]( 
Category 1 , Problem 1 ). The exact solution given there 

is: 


u = 0.5 exp 


-( Zn2)(^) 2 . 


In this example, Ax = 1 is chosen and At = 1 is based 
on CFL number = 1. With CFL number = 1, and other 
parameters e = 0 and a = 0, the CE/SE scheme yields a 
numerical result which is identical to the exact solution 
in the interior of the domain —20 ^ x ^ 451. Thus, 
performance of the outflow Type I NRBC at £ = 451 
can be easily validated. At t = 450, the Gaussian pulse 
is passing through the outflow boundary x = 451, where 
the Type I NRBC (12) with appropriate modification to 
1-D flow is imposed. The table in Fig. 8 lists the exact 
solution and the CE/SE result from x = 430 to x = 
451. The result is also plotted in Fig. 9. It is seen that 
they are completely identical (up to 10 decimal places) 
and there is absolutely no reflection, although the grid is 
rather crude with Ax = 1. 


3.2 2-D free shear layer instability and vortex 

roll-up 

The problem considered here is identical to the inviscid 
free shear layer instability problem considered in [14] 
(Fig. 10). The background mean flow consists of a fast 
stream (supersonic) in the upper half domain and a slow 


stream (subsonic) in the lower half domain. The two par- 
allel streams are connected by a continuously changing 
shear layer of the hyperbolic tangent profile. 

In the test, two computational domains are chosen. 
The first one is 0 ^ x 200 and —10 ^ y ^ 10, 
with a grid of 200x 100 uniform cells. The second one 
is 0 ^ x ^ 100 and — 10 ^ y ^ 10, with a grid of 
100 xl 00 uniform cells. Both cases have exactly the 
same grid cell sizes, time step size At = 0.15, and pa- 
rameters e = 0.2, a = 0. They are both run for 4000 
time steps when the spatial instability is fully developed. 
To ensure that the instability waves and vortex roll-up 
develop quickly, a large perturbation amplitude of 0.02 
at the most unstable frequency is chosen for the eigen- 
functions. At the outflow boundaries, the Type I outflow 
NRBC is used. Figure 1 1 demonstrates snapshots of the 
isobars and isopycnics in the two cases. These contours 
are observed to be almost identical to each other in their 
common domain portion. The contours in the short do- 
main seem as if they were a piece chopped off from the 
longer one. This shows that the outflow NRBC in this 
case is nearly perfect. 

3.3 Acoustic Pulse Propagation 

This problem is a typical subsonic wave propagation 
problem [13]. The computational domain in the x-y 
plane is a square with -100 ^ x ^ 100, and 0 ^ y ^ 
200. A uniform 201x 201 (triangulated) grid is used with 
Ax = Ay = 1 . Initially, a Gaussian acoustic pulse is lo- 
cated at the lower portion of the domain (x = 0, y = 25), 
with a mean flow of Mach number M = 0.5 in x direc- 
tion and a solid wall at the bottom. At the other three 
boundaries, Type II NRBC is used. By choosing a small 
amplitude factor S = 0.001, the Euler equations are prac- 
tically linearized. Fig. 12 shows the isobars at t = 100 
and the comparison between numerical and analytical re- 
sults for density along the line x = y for 0 ^ x 100. 
Although the wave propagation direction is oblique to the 
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200x20 domain, 200x100 grid 
(same contour levels) 

- 100x20 domain, 100x100 grid 


comparison of numerical results with different computational domains, 
showing effectiveness of outflow NRBC 




200x20 domain, 200x100 grid 
(same contour levels) 
100x20 domain, 100x100 grid 


Figure 1 1: Contours for long and short domains, showing effectiveness of the outflow NRBC. 



Figure 12: An acoustic pulse above a solid surface passing through the outflow boundary. 
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Figure 13: Sketch of the rectangular jet, aspect ratio 5, 
jet Mach number Mj= 1 . 6 , L = 1 6D, W = 1 6D and 
H = 5.6 D. 

outflow boundary, only tiny reflection is observed from 
Fig. 12. Another 2-D example with completely subsonic 
flow will be illustrated in §4. 

3.4 3-D rectangular jet flow 

Fig. 13 is the sketch of an underexpanded rectangular jet 
in 3-D space. The rectangular nozzle protrudes into the 
computational domain by l = 2D, D being the width of 
the jet. The unstructured mesh consists of about 1.7 mil- 
lion tetrahedral cells. At the inlet plane, ambient (station- 
ary) condition is specified. Jet flow at higher pressure is 
specified at the nozzle exit. All the rest boundaries are ei- 
ther Type I or Type II NRBC. Fig.s 1 4 shows snapshots of 
the isobars and v velocity contours on the cross sectional 
mid-planes after running 60,000 time steps. Fig. 15 
demonstrates the 3-D pressure iso-surfaces. No visible 
reflection is observed. 

3.5 Influence of the NRBC to the numerical 
accuracy 

When the outflow boundary is non-reflecting, the influ- 
ence from the boundary to the interior flow is small. 
Fig. 1 6 demonstrates the sound intensity level computed 
based on two domains with different lengths but same 
width for the 2-D Mach radiation problem in the 3rd 
CAA Workshop (Category 5) [ 16]. The short domain has 
a grid of 280 x 144 nodes. The only difference is that the 
longer domain has 30 more uniformly distributed nodes 
added in the x direction downstream. Fig. 16 shows the 
sound intensity levels ( square of r.m.s. p ' - pressure fluc- 
tuation ) along the line y = 10 at t = 400 (40000 time 
steps or 28 periods ) for both domains. It is observed 
that the maximum difference is about 2 x 10 -l °, far be- 
low the discretization error, thus is negligible. This case 
also demonstrates the relative drifting side effect of the 



v- contours on the mid plane (narrow side) 



v-contours on the mid-plane (wide side) 


Figure 14: v velocity contours on the mid-planes with 
mesh background. No visible reflection is observed. At 
the outflow boundary, flow is supersonic in the jet core 
and then becomes subsonic across the thick shear layer. 


NRBC 



Figure 15: Pressure iso-surfaces, no reflection is ob- 
served. 
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Figure 16: Comparison of sound intensity levels from 
computational domains of different lengths, showing that 
the outflow NRBC has negligible infuence to the interior. 

extrapolation Type I NRBC, but the error is acceptable. 

4 Numerical example for the absorbing 
boundary condition 

As shown in §2.4, for an absorbing boundary condi- 
tion such as the inlet, the prescribed inlet conditions are 
already good enough for an NRBC. Fig. 17 illustrates the 
instantaneous isobars for a cavity noise problem. 42, 000 
triangulated structured cells are used. The problem is 
a M=0.8 flow past a cavity of aspect ratio of 6.5. At 
the cavity walls, no slip boundary condition is imposed. 
Due to vortex shedding and acoustic feedback at the cav- 
ity edges, strong nonlinear acoustic waves are generated 
and propagate in both upstream and downstream direc- 
tions [15]. Fig. 18 is an enlargement of Fig. 17 around 
the inlet area. The details of the contours at the matched 
layer is revealed. It is observed that there is no spurious 
reflection and that the matched layer is about 4-5 cells 1 
thick. The matched layer in the Type II NRBC is some- 
what similar to that of the PML (perfectly matched layer) 
method [5,6] in that the difference diminishes quickly 
within this layer. But the layer arises automatically and 
there is no need to solve a new set of equations in the 
layer or to impose any conditions other than the pre- 
scribed inflow physical conditions. 

For NRBC at the top of the computational domain, 
a Type I NRBC may be modified from (12) by simply 
exchanging the axes x and y : Uq = U p ,(U y )g = 
(U X ) Q = 0, where as before, P and Q are respectively 
the boundary cell center and the ghost cell center at the 
top boundary. From Fig.s 17-18, even when the acoustic 
wave is oblique to the boundary, there is still no visi- 



Figure 17: Isobars snapshot for a cavity noise problem 
(Mach number M = 0.8), showing inflow NRBC and its 
absorbing property. No visible reflection is found at the 
top and the outflow boundary 

ble reflection. A Type I NRBC is applied to the outflow 
boundary, still, no visible reflection is found for the sub- 
sonic flow. 

5 Application of the buffer/sponge zone 

Generally, the NRBCs amount to very little reflection. 
However, there are situations that they may fail and dis- 
cernible reflections occur. A simple but effective rem- 
edy is to add a buffer/sponge zone between the core do- 
main and the boundary. In the buffer zone, the same 
governing equations are still used, except that the cell 
size in the buffer zone may grow rapidly (e.g. exponen- 
tially ) and create larger numerical damping. Typically, 
the number of cells in a buffer zone may vary from a 
few to 20. An example of a 2-D axisymmetric jet with 
Mach radiation from an externally stimulated shear layer 
is demonstrated. It is similar to the one described in the 
3rd CAA workshop benchmark problems (Category 5) 
[16]. The domain size is 33D x 19D with D being the 
jet nozzle diameter. 300 x 280 non-uniform rectangular 
cells are used before they are further triangulated. Ini- 
tially, a Mach number M = 2 jet exists. At the center 
of the nozzle exit plane, a source is imposed and per- 
turbs the jet flow with a small amplitude A = 0.001 at 
a Strouhal number St = 0.2. Mach radiation is then 
triggered and gains its strength along the stream. At the 
outflow boundary, particularly at the shear layer, due to 
the staggered type mesh and the strong Mach waves, the 
Type I NRBC fails and spurious reflection is generated 
and propagates upstream (Fig. 19). After a buffer zone 
of 10 cells is added at the domain outflow boundary, the 
spurious reflection disappears and a clean mach wave ra- 
diation is shown (Fig. 20). The size of the buffer zone 
cells grows exponentially at a rate of 207c along the x 
axis. 
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Figure 1 8: Details of the contours at the inflow boundary, 
showing the matched layer at the inlet and its spreading 
over the grid. 


Figure 19: Mach radiation from a M = 2 axisymmetric 
jet (without buffer zone), showing severe spurious nu- 
merical reflection. 



Figure 20: Mach radiation from a M = 2 axisymmetric 
jet (with buffer zone but not shown), showing a clean 
acoustic field. 

6 Concluding Remarks 

In the present paper, the hyperbolicity of the Euler 
equation system and the propagation of plane waves are 
revisited, and then combined to derive the continuity 
criterion of NRBC. The relation between flux balance 
across the boundary surface and the NRBC is estab- 
lished. Simple but effective C l continuity NRBCs are 
consequently developed. 

These NRBCs are simple and robust, as demonstrated 
in the one and multi-dimensional numerical examples 
based on the the recent CE/SE method. Numerical 
examples with other finite volume schemes, e.g., N.T. 
(Nessyahu and Tadmor) or upwind schemes can be found 
in [17) (§8). Generally, their performances are similar to 
those of the characteristics-based NRBCs. Limitations of 
the NRBCs are also discussed in §2.5. In particular, the 
Type I extrapolation NRBCs may cause solution drifting 
due to lack of information beyond the (outflow) bound- 
ary. A remedy is to incorporate the physical boundary 
conditions or to use a sponge zone. 

The diversity of various NRBCs in flow computations 
can never be overestimated. The purpose of the present 
paper is to show some guidelines in this direction and 
develop NRBCs that are simple but robust for practical 
computations. The compact updating procedure proves 
to work well with the NRBCs and provides C l accu- 
racy in surface flux evaluation. But it is definitely not 
the only way to achieve non-reflecting effect. Different 
schemes may have different treatments. Sometimes, a 
combination of the NRBC treatments may provide much 
improved results, such as the incorporation of the buffer 
zone. As a byproduct, we are now able to explain why 
the NRBCs with the recently developed CE/SE scheme 
[10,1 1,14,15] are robust. 

The restriction in §2.1 that the flow is continuous may 
be lifted, since for shock-capturing schemes, a discon- 
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tinuity may be considered as a continuous wave with 
steeper gradient. But if the wave propagation direction is 
oblique to the (artificial) boundary, the NRBC may per- 
form poorly due to reasons explained in §2.5. Chang et 
al [9], Huynh [17] presented 1-D examples ( i.e. wave 
propagation direction normal to the boundary) showing 
how a shock passes through an artificial boundary with- 
out causing visible reflection. 

At last, we comment on the influence of the errors of 
the NRBCs on the accuracy for domain interior. Gener- 
ally speaking, this involves the stability property of the 
scheme. If the time-marching explicit scheme satisfies 
the von Neumann stability criterion, an error introduced 
by the boundary condition will decay exponentially in 
time and space when it convects towards the domain in- 
terior. 
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Appendix I: Plane wave solutions 

The plane wave solutions are based on the Cauchy’s 

method of Fourier Integral (see Courant and Hilbertf 1 ], 

pp. 210-211) for linear homogeneous differential equa- 
tion. Consider Eq. (8): 


0V ,dV ~dV ~dV - 

~dt +A ~d^ + B ~d^ + c ~d^~ Q ' 
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Figure 21: Continuity criterion in 1-D flow. 


In the neighborhood of a point O (x 0 , to) at the artifi- 
cial boundary, (8) can be locally linearized by setting the 
jacobians A, B , C to their values at O. Assume a plane 
wave with the form V = We 16 , and substitute in (8). 
Here 0 = 0(x, t) = x • k — ujt is the phase of the simple 
wave, i = \f—\. 9 = const, stands for a characteristic 
surface or wave front [1]. (8) then becomes: 


,0V ,0V 6 0V ., v ~ 


where I is the 5 >6 identity matrix. For any given wave 
number vector k = ( k x ,k y ,k z ), from the hyperbolicity 
of (8), real cus exist such that K - uol = 0 (dispersion 
relation), where the matrix K — k x A - h k y B + k z C. 
Then the ‘amplitude’ V may be solved from 


0V ~0V 
dt + dx 




For more general waves other than the simple plane 
waves, as (8) is locally linearized in the neighborhood of 
the boundary point (9(x 0 , £ 0 )> they may be decomposed 
by Fourier integral with respect to wave number k and 
replaced by the superposition of plane waves. 


Appendix II: Discussions on the continuity 
criterion and the extrapolation technique 

Without loss of generality, consider a 1-D flow shown 
in Fig. 21. As a pure initial value Cauchy problem, 
V = V(ar, 0) at t = 0 is given. This problem is well- 
posed and there exists a unique soluiton V = V 0 (x, t) 
over the entire domain D\ — oo < x < oo, t ^ 0. 
If a boundary exists on the left inlet side, the problem 
becomes an initial-boundary-value problem, a boundary 
condition is required at the inlet boundary. 

Assume the artificial boundary locates at x = 0 and 
bisects the entire domain D into domain interior Di and 
domain exterior D e : D = D t + D e (see §2.1). In order 
that Vj and V e in §2.1 are respectively the unique well- 
posed solutions in subdomains Di and D e , an admissible 
common boundary condition they share at x = 0 is: 


V i (0,<) = v c (0 l t) = v o (0 1 «). 

Here, V 0 ((M) is the entire domain solution along the 
line x = 0; or equivalently, V o (0,£) may be obtained 


by the method of characteristics as sketched in Fig. 2 1 at 
point O at the boundary. 

However, for NRBC problems in reality, 
V 0 (x,t), x ^ 0; or V e (x,^) is totally missing 
(otherwise there is no need to investigate the NRBCs). 
In this situation, the extrapolation technique e.g. (12) 
is not an unreasonable choice for the NRBC continuity 
criterion. The initial-boundary value problem for V* is 
well-posed only when the flow is supersonic or is known 
to remain unchanged across the artificial boundary at 
x = 0 by a priori information. Otherwise, even though 
locally the extrapolation (12) leads to non-reflecting 
effect, globally the ‘solution’ will keep drifting away and 
deviate from the true solution. The following are some 
remedies for practical numerical flow computations: 

[ 1 ] use a sponge (buffer) zone with highly increased nu- 
merical damping to filter away the wave ingredients 
in the flow; when the Type I extrapolation NRBC 
is applied to the outer outflow boundary, the flow is 
already uniform at the boundary. 

[2] incorporate the extrapolation with other physical 
boundary conditions (e.g. back pressure etc.). 

Notice that the present discussions do not apply to the 
Type II NRBCs. 
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